/////////////////////////////////////////////////////////////////////
//
//
//
//
//
// Compile for linux (CUDA toolkit already installed)
// export PATH=/usr/local/cuda/bin:$PATH
// export LD_LIBRARY_PATH=/usr/local/cuda/lib64:$LD_LIBRARY_PATH
//
// nvcc -arch=sm_20 -O ninePixelFilter.cu
/////////////////////////////////////////////////////////////////////

#include <iostream>
#include <stdlib.h>
using namespace std;
#define defaultHeight 10
#define defaultWidth 10
#define threadsInBlock 144
#define maxThreads (defaultHeight*defaultWidth*3)

void CUDAninePixelFitlerWrapper(int *paddedPix1D, int* filteredPix1D);

int main(){
	int height = defaultHeight;
	int width = defaultWidth;
	int pixels[height][width];
	for (int i = 0; i < height; i++){
		for(int j = 0; j < width; j++){
			pixels[i][j] = rand()%256;
		}
	}
	for (int i = 0; i < height; i++){
		for(int j = 0; j < width; j++){
			cout << pixels[i][j] << "\t";
		}
		cout << endl;
	}

	int paddedPixels[height+2][width+2];
	// generate the padded picture by copying the outer row and columns.
	// first copy the original picture
	int i;
	int j;
	for (i = 1; i< height+1 ; i++){
		for (j = 1; j < width+1; j++){
			paddedPixels[i][j] = pixels[i-1][j-1];
		}
	}
	// next copy the rows and columns
	i = 0;
	for(j = 1; j < width+1; j++){
		paddedPixels[i][j] = pixels[i][j-1];
	}
	i = height+1;
	for(j = 1; j < width+1; j++){
		paddedPixels[i][j] = pixels[i-2][j-1];
	}
	j = 0;
	for(i = 1; i < height+1; i++){
		paddedPixels[i][j] = pixels[i-1][j];
	}
	j = width+1;
	for(i = 1; i < height+1; i++){
		paddedPixels[i][j] = pixels[i-1][j-2];
	}

	// Turn paddedPixels from a 2D array into a 1D array.
	int * paddedPix1D;
	int * filteredPix1D;
	paddedPix1D = (int *) malloc (sizeof(int) * (height+2) * (width+2));
	filteredPix1D = (int *) malloc (sizeof(int) * (height) * (width));

	for (int i = 0; i < height+2; i++){
		for (int j = 0; j < width+2; j++){
			paddedPix1D [i*(width+2) + j] = paddedPixels[i][j];
		}
	}

	CUDAninePixelFitlerWrapper(paddedPix1D, filteredPix1D);

	cout<< "filtered Pixels" << endl;
	for (int i = 0; i < height; i++){
		for(int j = 0; j < width; j++){
			cout << filteredPix1D[i*width+j] << "\t";
		}
		cout << endl;
	}

}

__global__ void kernelFilterNinePixels (int *dev_paddedPix1D, int *dev_filteredPix1D){

	int tid =threadIdx.x *blockDim.y + threadIdx.y;
	int blockID = blockIdx.x*gridDim.y + blockIdx.y;

	int overallTid = blockID*threadsInBlock + tid;
	int pixelTid = overallTid/3;
	// figure out which pixel each thread is working on
	// by decoding the (x,y) values
	int i = pixelTid / defaultWidth;
	int j = pixelTid % defaultWidth;
	int min, max, med, minMax, medMed, maxMin, finalMed;

	__shared__ int ninePixels[threadsInBlock*3];
	__shared__ int ninePixelsRowSorted[threadsInBlock*3];
	__shared__ int minMedMax[threadsInBlock];
	// EX: threads 0, 1, 2 of block 0 work with pixel (0,0)
	// thread 0 writes to ninePixels[0 1 2]
	// thread 1 writes to ninePixels[3 4 5]
	// thread 2 writes to ninePixels[6 7 8]

	// thread 0 fetches from paddedPix2D (0,0) (0,1) (0,2) = paddedPix1D[0 1 2]
	// thread 1 fetches from paddedPix2D (1,0) (1,1) (1,2)
	// = paddedPix1D[defaultWidth+0 defaultWidth+1 defaultWidth+2]
	// thread 2 fetches from paddedPix2d (2,0) (2,1) (2,2)
	// = paddedPix1D[2*defaultWidth+0 2*defaultWidth+1 2*defaultWidth+2]

	if (overallTid < maxThreads){
		if (tid % 3 == 0){
			ninePixels[tid*3] = dev_paddedPix1D[i*defaultWidth + j];
			ninePixels[tid*3+1] = dev_paddedPix1D[i*defaultWidth + j + 1];
			ninePixels[tid*3+2] = dev_paddedPix1D[i*defaultWidth + j + 2];

		}
		else if (tid % 3 == 1){
			ninePixels[tid*3] = dev_paddedPix1D[(i+1)*defaultWidth + j];
			ninePixels[tid*3+1] = dev_paddedPix1D[(i+1)*defaultWidth + j];
			ninePixels[tid*3+2] = dev_paddedPix1D[(i+1)*defaultWidth + j + 1];
		}
		else {
			ninePixels[tid*3] = dev_paddedPix1D[(i+2)*defaultWidth + j];
			ninePixels[tid*3+1] = dev_paddedPix1D[(i+2)*defaultWidth + j];
			ninePixels[tid*3+2] = dev_paddedPix1D[(i+2)*defaultWidth + j + 1];
		}
		__syncthreads();


		if (ninePixels[tid*3+0]> ninePixels[tid*3+1]){ // 0 > 1
			if (ninePixels[tid*3+0] > ninePixels[tid*3+2]){ // 0 > 2
				max = ninePixels[tid*3+0]; // max
				if (ninePixels[tid*3+1]>ninePixels[tid*3+2]){ // 0 > 1 > 2
					med = ninePixels[tid*3+1];
					min = ninePixels[tid*3+2];
				}
				else{ //  0 > 2 > 1
					med = ninePixels[tid*3+2];
					min = ninePixels[tid*3+1];
				}
			}
			else { // 2 > 0 > 1
				max = ninePixels[tid*3+2];
				med = ninePixels[tid*3+0];
				min = ninePixels[tid*3+1];
			}
		}

		else { // 1 > 0
			if (ninePixels[tid*3+1] > ninePixels[tid*3+2]){ // 1 > 2
				max = ninePixels[tid*3+1]; // max
				if (ninePixels[tid*3+0] > ninePixels[tid*3+2]){ // 1 > 0 > 2
					med = ninePixels[tid*3+0];
					min = ninePixels[tid*3+2];
				}
				else { // 1 > 2 > 0
					med = ninePixels[tid*3 + 2];
					min = ninePixels[tid*3 + 0];
				}
			}
			else { // 2 > 1 > 0
				max = ninePixels[tid*3+2];
				med = ninePixels[tid*3+1];
				min = ninePixels[tid*3+0];
			}
		}
		ninePixelsRowSorted[tid*3] = min;
		ninePixelsRowSorted[tid*3+1] = med;
		ninePixelsRowSorted[tid*3+2] = max;
		__syncthreads();

		// Find minMax of ninePixelsRowsSorted[tid*3+2]
		// minMedMax[0 1 2] will hold minMax, medMed, and maxMin

		// this is max of 3 Mins
		// thread 0 should find the max of ninePixelsRowSorted[0 3 6]
		// thread 3 should find the max of ninePixelsRowSorted[9 12 15], etc.
		if (tid % 3 == 0){
			if (ninePixelsRowSorted[tid*3] > ninePixelsRowSorted[(tid+1)*3]){ //0 > 1
				if (ninePixelsRowSorted[tid*3] > ninePixelsRowSorted[(tid+2)*3]) // 0 > 2
					minMedMax[tid] = ninePixelsRowSorted[tid*3];

				else //2>0>1
					minMedMax[tid] = ninePixelsRowSorted[(tid+2)*3];

			}
			else { //1 > 0
				if (ninePixelsRowSorted[(tid+1)*3] > ninePixelsRowSorted[(tid+2)*3]) // 1 > 2
					minMedMax[tid] = ninePixelsRowSorted[(tid+1)*3];
				else // 2 > 1 > 0
					minMedMax[tid] = ninePixelsRowSorted[(tid+2)*3];
			}
		}

		// find the med of 3 medians
		// thread 1 should find the med of ninePixelsRowSorted[1 4 7]
		// thread 4 should find the med of ninePixelsRowSorted[10 13 16]
		else if (tid % 3 == 1){
			if (ninePixelsRowSorted[tid*3+1] > ninePixelsRowSorted[(tid+1)*3 + 1]){ // 0 > 1
				if (ninePixelsRowSorted[tid*3 + 1] > ninePixelsRowSorted[(tid+2)*3 + 1]){ // 0 > 2
					if (ninePixelsRowSorted[(tid+1)*3 + 1] > ninePixelsRowSorted[(tid+2)*3 + 1]) // 0 > 1 > 2
						minMedMax[tid] = ninePixelsRowSorted[(tid+1)*3 + 1] ;
					else // 0 > 2 > 1
						minMedMax[tid] = ninePixelsRowSorted[(tid+2)*3 + 1] ;
				}
				else // 2 > 0 > 1
					minMedMax[tid] = ninePixelsRowSorted[tid*3 + 1] ;
			}
			else{ //1 > 0
				if (ninePixelsRowSorted[(tid+1)*3+1] > ninePixelsRowSorted[(tid+2)*3+1]){ // 1> 2
					if (ninePixelsRowSorted[(tid+2)*3+1] > ninePixelsRowSorted[tid*3+1]) // 1 > 2 > 0
						minMedMax[tid] = ninePixelsRowSorted[(tid+2)*3 + 1] ;
					else // 1 > 0 > 2
						minMedMax[tid] = ninePixelsRowSorted[tid*3 + 1];
				}
				else //  2 > 1 > 0
					minMedMax[tid] = ninePixelsRowSorted[(tid+1)*3 + 1] ;

			}
		}

		// find the min of 3 maximums
		// thread 2 should find the min of ninePixelsRowSorted[2 5 8]
		// thread 5 should find the med of ninePixelsRowSorted[11 14 17]

		else{ // (tid %3 == 2)
			if (ninePixelsRowSorted[tid*3] < ninePixelsRowSorted[(tid+1)*3]){ //0 < 1
				if (ninePixelsRowSorted[tid*3] < ninePixelsRowSorted[(tid+2)*3]) // 0 < 2
					minMedMax[tid] = ninePixelsRowSorted[tid*3];

				else // 2 < 0 < 1
					minMedMax[tid] = ninePixelsRowSorted[(tid+2)*3];

			}
			else { //1 < 0
				if (ninePixelsRowSorted[(tid+1)*3] < ninePixelsRowSorted[(tid+2)*3]) // 1 < 2
					minMedMax[tid] = ninePixelsRowSorted[(tid+1)*3];
				else // 2 < 1 < 0
					minMedMax[tid] = ninePixelsRowSorted[(tid+2)*3];
			}
		}

		__syncthreads();

		// Now find the final median value
		// use 1 in 3 threads to do this (using threads 0,3,6,9, etc)
		// thread 0 finds median of minMedMax[0 1 2]
		// thread 3 finds the median of minMedMax[3 4 5]
		// thread 0 needs to write to dev_filteredPix1D[pixelTid]
		// note that threads 0,1,2 of block 0 have the same pixelTid
		// but the pixelTid of thread 0 of block 1 = pixelTid block 0 thread 0 + 48.


		if (tid %3 == 0){
			if (minMedMax[tid] >minMedMax[tid+1]){ // 0 > 1
				if (minMedMax[tid] > minMedMax[tid+2]){
					if (minMedMax[tid+1] > minMedMax[tid+2]) // 0 > 1 > 2
						dev_filteredPix1D[pixelTid]= minMedMax[tid+1] ;
					else // 0 > 2 > 1
						dev_filteredPix1D[pixelTid]= minMedMax[tid+2] ;
				}
				else{ // 2 > 0 > 1
					dev_filteredPix1D[pixelTid] = minMedMax[tid];
				}
			}

			else{ //1 > 0
				if (minMedMax[tid+1]  > minMedMax[tid+2] ){ // 1> 2
					if (minMedMax[tid+2] > minMedMax[tid]) // 1 > 2 > 0
						dev_filteredPix1D[pixelTid]= minMedMax[tid+2] ;
					else // 1 > 0 > 2
						dev_filteredPix1D[pixelTid]= minMedMax[tid] ;
				}
				else{ //  2 > 1 > 0
					dev_filteredPix1D[pixelTid]=minMedMax[tid+1];
				}
			}
		}

	}// end of if(overallTid < maxThreads)

} // end of kernel








////////////////////////////////////////////////////////////////////
//
// This wrapper calls sets up calls to the CUDA enabled device.
//
// The wrapper uses global memory and 1 dimensional arrays.
// The CUDA algorithm implemented is a median filter, with 3x3 pixels
// This wrapper will NOT work for different sized median filters.
//
///////////////////////////////////////////////////////////////////

void CUDAninePixelFitlerWrapper(int *paddedPix1D, int* filteredPix1D){
	int height = defaultHeight;
	int width = defaultWidth;
	int numPixels = height*width;
	int *dev_paddedPix1D, *dev_filteredPix1D;

	// allocate memory
	cudaMalloc( (void **)&dev_paddedPix1D, sizeof(int) * (height+2) * (width+2) );
    cudaMalloc( (void **)&dev_filteredPix1D, sizeof(int) * (height) * (width) );

    // copy padded pixels from host to device
    cudaMemcpy(dev_paddedPix1D, paddedPix1D, sizeof(int)*(height+2)*(width+2), cudaMemcpyDeviceToHost );

    // figure out the dimensions of the CUDA kernel launch
    // The scheme is to launch 3 threads per pixel
    // with 144 threads = 48 pixels per block
    // figure out the number of blocks to launch
    int numberofBlocks = (numPixels%48 == 0) ? numPixels/48 : numPixels/48+1;

    // figure out the block configuration
    int numBlocksX, numBlocksY;
    numBlocksY = numberofBlocks % 16 + 1;
    numBlocksX = numberofBlocks / 16 + 1;

    dim3 Grid(numBlocksX, numBlocksY);
    dim3 Block(16, 9);	// set at 144 threads

    kernelFilterNinePixels <<<Grid, Block>>> (dev_paddedPix1D, dev_filteredPix1D);

    // Copy filtered pixels from device back to host.
    cudaMemcpy(filteredPix1D, dev_filteredPix1D, sizeof(int)*height*width, cudaMemcpyDeviceToHost );
    return;
}
